home *** CD-ROM | disk | FTP | other *** search
/ FishMarket 1.0 / FishMarket v1.0.iso / fishies / 376-400 / disk_386 / xlispstat / src2.lzh / XLisp-Stat / splines.c < prev    next >
C/C++ Source or Header  |  1990-10-04  |  4KB  |  122 lines
  1. /* natural cubic spline interpolation based on Numerical Recipes in C */
  2. /* XLISP-STAT 2.1 Copyright (c) 1990, by Luke Tierney                  */
  3. /* Additions to Xlisp 2.1, Copyright (c) 1989 by David Michael Betz    */
  4. /* You may give out copies of this software; for conditions see the    */
  5. /* file COPYING included with this distribution.                       */
  6.  
  7. #include "xlisp.h"
  8. #include "osdef.h"
  9. #ifdef ANSI
  10. #include "xlsproto.h"
  11. #else
  12. #include "xlsfun.h"
  13. #endif ANSI
  14.  
  15. #ifdef ANSI
  16. void find_spline_derivs(double *,double *,int,double *,double *),
  17.      spline_interp(double *,double *,double *,int,double,double *);
  18. #else find_spline_derivs(),
  19.      spline_interp();
  20. #endif ANSI
  21.  
  22. /* calculate second derivatives; assumes strictly increasing x values */
  23. static void find_spline_derivs(x, y, n, y2, u)
  24.      double *x, *y, *y2, *u;
  25.      int n;
  26. {
  27.   int i, k;
  28.   double p, sig;
  29.  
  30.   y2[0] = u[0] = 0.0;  /* lower boundary condition for natural spline */
  31.   
  32.   /* decomposition loop for the tridiagonal algorithm */
  33.   for (i = 1; i < n - 1; i++) {
  34.     y2[i] = u[i] = 0.0; /* set in case a zero test is failed */
  35.     if (x[i - 1] < x[i] && x[i] < x[i + 1]) {
  36.       sig = (x[i] - x[i - 1]) / (x[i + 1] - x[i - 1]);
  37.       p = sig * y2[i - 1] + 2.0;
  38.       if (p != 0.0) {
  39.         y2[i] = (sig - 1.0) / p;
  40.         u[i] = (y[i + 1] - y[i]) / (x[i + 1] - x[i])
  41.              - (y[i] - y[i - 1]) / (x[i] - x[i - 1]);
  42.         u[i] = (6.0 * u[i] / (x[i + 1] - x[i - 1]) - sig * u[i - 1]) / p;
  43.       }
  44.     }
  45.   }
  46.  
  47.   /* upper boundary condition for natural spline */
  48.   y2[n - 1] = 0.0;
  49.  
  50.   /* backsubstitution loop of the tridiagonal algorithm */
  51.   for (k = n - 2; k >= 0; k--)
  52.     y2[k] = y2[k] * y2[k + 1] + u[k];
  53. }
  54.  
  55. /* interpolate or extrapolate value at x using results of find_spline_derivs */
  56. static void spline_interp(xa, ya, y2a, n, x, y)
  57.      double *xa, *ya, *y2a, x, *y;
  58.      int n;
  59. {
  60.   int klo, khi, k;
  61.   double h, b, a;
  62.  
  63.   if (x <= xa[0]) {
  64.     h = xa[1] - xa[0];
  65.     b = (h > 0.0) ? (ya[1] - ya[0]) / h - h * y2a[1] / 6.0 : 0.0;
  66.     *y = ya[0] + b * (x - xa[0]);
  67.   }
  68.   else if (x >= xa[n - 1]) {
  69.     h = xa[n - 1] - xa[n - 2];
  70.     b = (h > 0.0) ? (ya[n - 1] - ya[n - 2]) / h + h * y2a[n - 2] / 6.0 : 0.0;
  71.     *y = ya[n - 1] + b * (x - xa[n - 1]);
  72.   }
  73.   else {
  74.     /* try a linear interpolation for equally spaced x values */
  75.     k = (n - 1) * (x - xa[0]) / (xa[n - 1] - xa[0]);
  76.  
  77.     /* make sure the range is right */
  78.     if (k < 0) k = 0;
  79.     if (k > n - 2) k = n - 2;
  80.  
  81.     /* bisect if necessary until the bracketing interval is found */
  82.     klo = (x >= xa[k]) ? k : 0;
  83.     khi = (x < xa[k + 1]) ? k + 1 : n - 1;
  84.     while (khi - klo > 1) {
  85.       k = (khi + klo) / 2;
  86.       if (xa[k] > x) khi = k;
  87.       else klo = k;
  88.     }
  89.  
  90.     /* interpolate */
  91.     h = xa[khi] - xa[klo];
  92.     if (h > 0.0) {
  93.       a = (xa[khi] - x) / h;
  94.       b = (x - xa[klo]) / h;
  95.       *y = a * ya[klo] + b * ya[khi]
  96.          + ((a * a * a - a) * y2a[klo] + (b * b * b - b) * y2a[khi]) * (h * h) / 6.0;
  97.     }
  98.     else *y = (ya[klo] + ya[khi]) / 2.0; /* should not be needed */
  99.   }
  100. }
  101.  
  102. int fit_spline(n, x, y, ns, xs, ys, work)
  103.      int n, ns;
  104.      RVector x,y,xs,ys,work;/* double *x, *y, *xs, *ys, *work; changed JKL */
  105. {
  106.   int i;
  107.   double *y2, *u;
  108.  
  109.   y2 = work; u = work + n;
  110.  
  111.   if (n < 2 || ns < 1) return (1); /* signal an error */
  112.   for (i = 1; i < n; i++)
  113.     if (x[i - 1] >= x[i]) return(1); /* signal an error */
  114.   
  115.   find_spline_derivs(x, y, n, y2, u);
  116.   
  117.   for (i = 0; i < ns; i++)
  118.     spline_interp(x, y, y2, n, xs[i], &ys[i]);
  119.  
  120.   return(0);
  121. }
  122.